Fire Parameterization

Ignitions

library(plyr)
library(pscl)

## Warning: package 'pscl' was built under R version 4.0.3

## Classes and Methods for R developed in the
## Political Science Computational Laboratory
## Department of Political Science
## Stanford University
## Simon Jackman
## hurdle and zeroinfl functions by Achim Zeileis

library(sf)

## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1

library(sp)
library(raster)

w_dir<-"C:/Users/zacha/Desktop/Sapps_DM_paper/"
CNF<-st_read("C:/Users/zacha/Desktop/Sapps_DM_paper/Georgia.shp")

## Reading layer `Georgia' from data source `C:\Users\zacha\Desktop\Sapps_DM_paper\Georgia.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 14 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 79536.62 ymin: 3358689 xmax: 520501.5 ymax: 3881045
## projected CRS:  NAD83 / UTM zone 17N

l_fire_dat <- read.csv(paste(w_dir,"Inputs/FiresInAppsLightning2.csv", sep=""))
colnames(l_fire_dat)

## [1] "X"           "FIRE_YEAR"   "DISCOVER_1"  "STAT_CAU_1"  "CONT_DATE"  
## [6] "FIRE_SIZE"   "LATITUDE"    "LONGITUDE"   "Coordinates"

#UTM11<-Plot_locations[Plot_locations$UTM.Zone==('11S'),]
xy <-l_fire_dat[,c("LONGITUDE","LATITUDE")]

Spati <- SpatialPointsDataFrame(coords = xy, data = l_fire_dat,
                               proj4string = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs "))%>%
  spTransform(projection(CNF))%>%
  as("sf")
Int<-st_intersection(Spati,CNF)

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

#as(Spati,)
plot(Spati)

l_fire_dat <- Int
#head(l_fire_dat)
l_fire_dat<-l_fire_dat%>%
  subset(FIRE_SIZE>15.44)
##This moves us from 676 to 155

l_fire_days <- as.data.frame(cbind(l_fire_dat$FIRE_YEAR, l_fire_dat$DISCOVER_1)) #Extracting year and julian day
colnames(l_fire_days) <- c("YEAR", "J_DAY")
l_fire_days_sort <- l_fire_days[order(l_fire_days[,1]),] #sorting by year

library(RColorBrewer)

red<-RColorBrewer::brewer.pal(9,'YlOrRd')
#plot no of fires/yr
l_fire_days_sort_count<- cbind(l_fire_days_sort, rep(1, nrow(l_fire_days_sort)))
l_fires_count <- cbind.data.frame(unique(l_fire_days_sort_count$YEAR) , 
                                 tapply(l_fire_days_sort_count$`rep(1, nrow(l_fire_days_sort))`,l_fire_days_sort_count$YEAR, sum))
colnames(l_fires_count) <- c("YEAR", "COUNT")

barplot(l_fires_count$COUNT, main ="No of ign/yr Lightning Observed",col=red,names.arg=l_fires_count$YEAR)

h_fire_dat <- read.csv(paste(w_dir, "Inputs/FiresInAppsHuman2.csv", sep=""))

xy <-h_fire_dat [,c("LONGITUDE","LATITUDE")]

Spati <- SpatialPointsDataFrame(coords = xy, data = h_fire_dat ,
                               proj4string = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs "))%>%
  spTransform(projection(CNF))%>%
  as("sf")
Int<-st_intersection(Spati,CNF)

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

h_fire_dat <-Int
h_fire_dat<-h_fire_dat%>%
  subset(FIRE_SIZE>15.44)
##Goes from 24789 to 1531
#head(h.fire.dat)
h_fire_days <- as.data.frame(cbind(h_fire_dat$FIRE_YEAR, h_fire_dat$DISCOVER_1)) #Extracting year and julian day
colnames(h_fire_days) <- c("YEAR", "J_DAY")
h_fire_days_sort <- h_fire_days[order(h_fire_days[,1]),] #sorting by year
#plot no of fires/yr
h_fire_days_sort_count<- cbind(h_fire_days_sort, rep(1, nrow(h_fire_days_sort)))
h_fires_count <- cbind.data.frame(unique(h_fire_days_sort_count$YEAR) , 
                                  tapply(h_fire_days_sort_count$`rep(1, nrow(h_fire_days_sort))`,h_fire_days_sort_count$YEAR, sum))
colnames(h_fires_count) <- c("YEAR", "COUNT")
barplot(h_fires_count$COUNT, main ="No of ign/yr Human Accidental Observed",col=red,names.arg=h_fires_count$YEAR)

ign_types <- c("Lightning", "HumanAccidental")

fire_days_list <- list(l_fire_days_sort, h_fire_days_sort) ##organizing all ignition types into a list


##Import daily historic FWI data
FWI_dat <- read.csv(paste0(w_dir,"Inputs/GA_Climate.csv"))

FWI_dat<-with(FWI_dat,aggregate(FWI,by=list(Year=Year,Timestep=Timestep),FUN=mean))


#l_fires_count$
colnames(FWI_dat) <-c("Year","Timestep","FWI")
#colnames(FWI_dat)[3]<-"FWI"
FWI_dat$ID <- paste(FWI_dat$Year, "_", FWI_dat$Timestep, sep="") #creating date identifier out of date and julian day
FWI_dat$Dates<-as.Date(FWI_dat$Timestep, origin=as.Date(paste0(FWI_dat$Year,"-01-01")))
#l_fire_count$Date<-as.Date(l_fire_count$J_DAY, origin=as.Date(paste0(l_fire_count$YEAR,"-01-01")))
fire_days_count <- ddply(l_fire_days_sort, .(l_fire_days_sort$YEAR, l_fire_days_sort$J_DAY), nrow) #
fire_days_count$Date<-as.Date(fire_days_count$`l_fire_days_sort$J_DAY`,origin=as.Date(paste0(fire_days_count$`l_fire_days_sort$YEAR`,"-01-01")))

#plot(FWI_dat$Dates,FWI_dat$FWI,type="l")
#points(fire_days_count$Date,fire_days_count$V1,col="red")

igns_list <- list()

for (i in 1:length(ign_types[1:2])){#THIS DOESN'T INCLUDE RX BURNS BUT THATS CAUSE WE ARE PROVIDING THOSE TO SCRAPPLE DIRECT
  ign_type_select <- fire_days_list[[i]] ##selecting each ignition type individually
  fire_days_count <- ddply(ign_type_select, .(ign_type_select[,1], ign_type_select[,2]), nrow) #finds duplicate rows in fire data
  fire_days_count$Dates<-as.Date(fire_days_count[,2],origin=as.Date(paste0(fire_days_count[,1],"-01-01")))
  ##Merging dataframes by year and julian day
  fire_days_short <- subset(fire_days_count, fire_days_count[,1] > 1993 & fire_days_count[,1] < 2016) ##restricting fire records to climate file years
  FWI_short <- subset(FWI_dat, FWI_dat$Year > 1993 & FWI_dat$Year < 2016) #restricting climate data to fire history records
  #plot(FWI_short$Dates,FWI_short$FWI,type="l",main =ign_types[i],ylab="Fire Weather Index")
  #points(fire_days_count$Date,fire_days_count$V1,col=adjustcolor("red",alpha.f=.3),pch=19,cex=2.0)
 #fire_days_short$
  FWI_fire_merge <- join(FWI_short, fire_days_short, by="Dates",type="left") ##Merging based on unique date id
  FWI_fire_number <- FWI_fire_merge[,c(1,3,8)] #pulling out FWI and number of fires
  FWI_fire_number[is.na(FWI_fire_number)] <- 0 #converting NAs to 0, meaning 0 fires
  par(mfrow=c(1,1))
  plot(FWI_fire_number[,2], FWI_fire_number[,3], main =ign_types[i], xlab = "dailyFWI", ylab = "noFires",cex=2.0) #plotting FWI against no of fires just to look at pattern
  igns_list[[i]] <- FWI_fire_number 
}

Lightning<-as.data.frame(igns_list[[1]])
colnames(Lightning)[3]<-"No_FIRES"

summary(lm(Lightning$No_FIRES~Lightning$FWI))

## 
## Call:
## lm(formula = Lightning$No_FIRES ~ Lightning$FWI)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.01875 -0.00641 -0.00438 -0.00181  2.99408 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -0.0030925  0.0020394  -1.516    0.129    
## Lightning$FWI  0.0005824  0.0001476   3.945 8.05e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0755 on 8028 degrees of freedom
## Multiple R-squared:  0.001935,   Adjusted R-squared:  0.00181 
## F-statistic: 15.56 on 1 and 8028 DF,  p-value: 8.051e-05

#Lightning<-Lightning[-1936,]
#Lightning$No_FIRES
###Lightning
zeroinf_mod <- zeroinfl(as.numeric(No_FIRES)~as.numeric(FWI),data=Lightning, dist="poisson")
summary(zeroinf_mod)

## 
## Call:
## zeroinfl(formula = as.numeric(No_FIRES) ~ as.numeric(FWI), data = Lightning, 
##     dist = "poisson")
## 
## Pearson residuals:
##      Min       1Q   Median       3Q      Max 
## -0.12623 -0.06901 -0.04564 -0.02442 37.70623 
## 
## Count model coefficients (poisson with log link):
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      2.02999    0.86284   2.353   0.0186 *  
## as.numeric(FWI) -0.21642    0.05499  -3.936  8.3e-05 ***
## 
## Zero-inflation model coefficients (binomial with logit link):
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     10.24538    1.28792   7.955 1.79e-15 ***
## as.numeric(FWI) -0.40606    0.09149  -4.438 9.06e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Number of iterations in BFGS optimization: 14 
## Log-likelihood: -191.7 on 4 Df

zero<-predict(zeroinf_mod,data.frame(FWI=Lightning$FWI), type = "zero")
#print(zero)
#n<-data.frame(FWI=Lightning$FWI)
lambda <- predict(zeroinf_mod,data.frame(FWI=Lightning$FWI), type = "count")
Simulation<-ifelse(rbinom(8030,size = 1, prob = zero) > 0, 0, rpois(8030, lambda = lambda))

ReturnDf<-cbind(Simulation,Lightning)
Sample<-aggregate(ReturnDf$Simulation,by=list(year=ReturnDf$Year),FUN=sum)
par(mfrow=c(1,2))
barplot(Sample$x, main ="No of ign/yr Lightning Simulated",col='black', ylim=c(0,10),names.arg=Sample$year)
#abline(v=15.7)
barplot(l_fires_count$COUNT, main ="No of ign/yr Lightning Observed",col=red,ylim=c(0,10),names.arg=l_fires_count$YEAR)

#abline(v=15.7)

FireWeathers<-seq(0,44,.1)
par(mfrow=c(1,2))
zeros<-predict(zeroinf_mod,data.frame(FWI=FireWeathers), type = "zero")
plot(FireWeathers,zeros,ylab="Liklihood of a zero")
lambda <- predict(zeroinf_mod,data.frame(FWI=FireWeathers), type = "count")
plot(FireWeathers,lambda,ylab="Number of Igntions")

par(mfrow=c(1,1))
FireWeathers<-rep(seq(0,50,.1),100)
lambda <- predict(zeroinf_mod,data.frame(FWI=FireWeathers), type = "count")
zeros<-predict(zeroinf_mod,data.frame(FWI=FireWeathers), type = "zero")
Simulation<-ifelse(rbinom(50100,size = 1, prob = zeros) > 0, 0, rpois(50100, lambda = lambda))
plot(FireWeathers,Simulation,col=adjustcolor("black",alpha.f = .01),pch=19,ylab=("Ignitions in a day"),xlab="FireWeather Index",
     main="Fire Weather Index Scale 100 Simulations")

Accidental<-as.data.frame(igns_list[[2]])
colnames(Accidental)[3]<-"No_FIRES"
summary(lm(Accidental$No_FIRES~Accidental$FWI))

## 
## Call:
## lm(formula = Accidental$No_FIRES ~ Accidental$FWI)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.0419 -0.0321 -0.0304 -0.0283  3.9733 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    0.0245757  0.0057240   4.293 1.78e-05 ***
## Accidental$FWI 0.0004618  0.0004143   1.115    0.265    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2119 on 8028 degrees of freedom
## Multiple R-squared:  0.0001547,  Adjusted R-squared:  3.02e-05 
## F-statistic: 1.242 on 1 and 8028 DF,  p-value: 0.265

###Lightning
Acc_zeroinf_mod <- zeroinfl(No_FIRES~FWI,data=Accidental, dist="poisson")
summary(Acc_zeroinf_mod)

## 
## Call:
## zeroinfl(formula = No_FIRES ~ FWI, data = Accidental, dist = "poisson")
## 
## Pearson residuals:
##     Min      1Q  Median      3Q     Max 
## -0.1636 -0.1481 -0.1455 -0.1421 20.5314 
## 
## Count model coefficients (poisson with log link):
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept) -0.94748    0.31055  -3.051  0.00228 **
## FWI          0.01392    0.02044   0.681  0.49577   
## 
## Zero-inflation model coefficients (binomial with logit link):
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  2.6605450  0.3217303   8.269   <2e-16 ***
## FWI         -0.0002711  0.0214416  -0.013     0.99    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Number of iterations in BFGS optimization: 13 
## Log-likelihood: -1041 on 4 Df

#Tst<-predict(zeroinf_mod,type="zero")
zero<-predict(Acc_zeroinf_mod,data.frame(FWI=Accidental$FWI), type = "zero")
#print(zero)
#n<-data.frame(FWI=Lightning$FWI)
lambda <- predict(Acc_zeroinf_mod,data.frame(FWI=Accidental$FWI), type = "count")
Simulation<-ifelse(rbinom(8030,size = 1, prob = zero) > 0, 0, rpois(8030, lambda = lambda))

ReturnDf<-cbind(Simulation,Lightning)
Sample<-aggregate(ReturnDf$Simulation,by=list(year=ReturnDf$Year),FUN=sum)
par(mfrow=c(1,2))
barplot(Sample$x, main ="No of ign/yr Accidental_Simulated",col='black', ylim=c(0,50),names.arg=Sample$year)

barplot(h_fires_count$COUNT, main ="No of ign/yr Human Accidental Observed",col=red,ylim=c(0,50),names.arg=h_fires_count$YEAR)

FireWeathers<-seq(0,50,.1)
par(mfrow=c(1,2))
zeros<-predict(Acc_zeroinf_mod,data.frame(FWI=FireWeathers), type = "zero")
plot(FireWeathers,zeros,ylab="Liklihood of a zero",ylim=c(0,1))
lambda <- predict(Acc_zeroinf_mod,data.frame(FWI=FireWeathers), type = "count")
plot(FireWeathers,lambda,ylab="Number of Igntions")

par(mfrow=c(1,1))
FireWeathers<-rep(seq(0,50,.1),100)
lambda <- predict(Acc_zeroinf_mod,data.frame(FWI=FireWeathers), type = "count")
zeros<-predict(Acc_zeroinf_mod,data.frame(FWI=FireWeathers), type = "zero")

Simulation<-ifelse(rbinom(50100,size = 1, prob = zeros) > 0, 0, rpois(50100, lambda = lambda))
plot(FireWeathers,Simulation,col=adjustcolor("black",alpha.f = .01),pch=19,ylab=("Ignitions in a day"),xlab="FireWeather Index",
     main="Fire Weather Index Scale 100 Simulations")

library(geosphere)

## Warning: package 'geosphere' was built under R version 4.0.3

library(maptools)

## Checking rgeos availability: TRUE

library(raster)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:raster':
## 
##     intersect, select, union

## The following objects are masked from 'package:plyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(rgdal)

## rgdal: version: 1.5-16, (SVN revision 1050)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.0.4, released 2020/01/28
## Path to GDAL shared files: C:/Users/zacha/Documents/R/win-library/4.0/rgdal/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 6.3.1, February 10th, 2020, [PJ_VERSION: 631]
## Path to PROJ shared files: C:/Users/zacha/Documents/R/win-library/4.0/rgdal/proj
## Linking to sp version:1.4-2
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.

library(rgeos)

## rgeos version: 0.5-5, (SVN revision 640)
##  GEOS runtime version: 3.8.0-CAPI-1.13.1 
##  Linking to sp version: 1.4-2 
##  Polygon checking: TRUE

library(geosphere)
library(nplr)

## Warning: package 'nplr' was built under R version 4.0.3

## 
## Attaching package: 'nplr'

## The following object is masked from 'package:raster':
## 
##     overlay

library(pscl)
library(ROCR)

## Warning: package 'ROCR' was built under R version 4.0.3

library(likelihood)

## Warning: package 'likelihood' was built under R version 4.0.3

## Loading required package: nlme

## 
## Attaching package: 'nlme'

## The following object is masked from 'package:dplyr':
## 
##     collapse

## The following object is masked from 'package:raster':
## 
##     getData

library(hexbin)
library(RColorBrewer)

w_dir <- "C:/Users/zacha/Desktop/Truncated_Scrapple/Inputs/"
out_dir <-  "C:/Users/zacha/Desktop/Truncated_Scrapple/Outputs/"

Here we are looking to fit a binomial model as to how probability of spread is related to the FWI, wind, and fuels.

For this we have wind data and Fire weather index data located to climate regions in LANDIS-II This is Gridmet data from the University of Idaho http://www.climatologylab.org/gridmet.html meta

The ecoregion map is used to locate the windspeed and the fwi data for each time step.

fwi_dat <- read.csv(paste0(w_dir,"11_Eco_Climate.csv"),stringsAsFactors = FALSE)
wind_map<- raster::raster(paste(w_dir,"11_Ecoregions.tif",sep=""))

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
##  but +towgs84= values preserved

wsv_dat <- read.csv(paste(w_dir,"Wind_1028.csv",sep=""),stringsAsFactors = FALSE)

colnames(wsv_dat) <- c("date", "eco1", 'eco2', 'eco3', 'eco4','eco5','eco6','eco8','eco9','eco10','eco11')

plot(wind_map)

We also have a fuel map representing landscape level fuels for the study area as well as an uphill azimuth map and a slope map (used to calculate the direction of spread)

ExampleRaster<-raster(paste(w_dir,"11_Ecoregions.tif",sep=""))

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
##  but +towgs84= values preserved

projy<-crs(ExampleRaster)

### This is the landis-II fuel imput for timestep 10
fccs_raster<-raster(paste(w_dir,"Fuel_129.IMG",sep=""))
#plot(fccs_raster)
crs(fccs_raster)<-projy
extent(fccs_raster)<-extent(ExampleRaster)
plot(fccs_raster)

Given the sites that are turned off on the LANDIS-II map, we interpolated between simulated fuel values to get an even coverage of fuels.

r=fccs_raster
r[r==0] <- NA

fill.na <- function(x, i=5) {
  if( is.na(x)[i] ) {
    return( round(mean(x, na.rm=TRUE),0) )
  } else {
    return( round(x[i],0) )
  }
}  

r2 <- focal(r, w = matrix(1,3,3), fun = fill.na, 
            pad = TRUE, na.rm = FALSE )
plot(r2)
r3<-focal(r2, w = matrix(1,3,3), fun = fill.na, 
            pad = TRUE, na.rm = FALSE )

plot(r3)

r4<-focal(r3, w = matrix(1,3,3), fun = fill.na, 
            pad = TRUE, na.rm = FALSE )

plot(r4)
fuel_map<-r4

plot(fuel_map)
writeRaster(fuel_map,paste0(w_dir,"Interpolated_FuelMap.tif"))

fuel_map<-raster(paste0(w_dir,"Interpolated_FuelMap.tif"))

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
##  but +towgs84= values preserved

uphill_azi_map<- raster(paste(w_dir,"SAPPS_aspect_inverted_atModelRes.tif",sep=""))

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
##  but +towgs84= values preserved

#uphill_azi_map <- projectRaster(slope_azi_map, to = fuel_map, method='bilinear')
slope_map <-  raster(paste(w_dir,"SAPPS_Slope_ChrisArcmapprocess_modelres.tif",sep=""))

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
##  but +towgs84= values preserved

#
climate_stack <- stack(wind_map, fuel_map, uphill_azi_map, slope_map)
plot(climate_stack)

#Empty objects and column names and stuff used in loops
newfire<-readOGR(paste(w_dir,"fire_geomac_loc.shp",sep=""))

## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\zacha\Desktop\Truncated_Scrapple\Inputs\fire_geomac_loc.shp", layer: "fire_geomac_loc"
## with 591 features
## It has 64 fields
## Integer64 fields read as strings:  OBJECTID fireyear

## Putting into sets based on differnces in meta-data (columns not matched)
newfire1<-newfire[newfire$incidentna!=0,]
newfire2<-newfire[newfire$incidentna==0,]
newfire_df<-data.frame(newfire)
perimeter_map<-newfire
### Here we are cleaning the data
perimeter_map$TrueDate<-as.Date(perimeter_map$perimeterd)
perimeter_map$TrueDate[is.na(perimeter_map$TrueDate)]<-as.Date(perimeter_map$date_[is.na(perimeter_map$TrueDate)])
### Subseting to our time period of climate coverage (< 2017)
perimeter_map<-perimeter_map[perimeter_map$TrueDate<"2017-01-01",]
perimeter_map$incidentna[perimeter_map$incidentna=='0']<-perimeter_map$fire_name[perimeter_map$incidentna=='0']
perimeter_map$incidentna[perimeter_map$incidentna==0]<-perimeter_map$fire_name[perimeter_map$incidentna==0]

DFP<-as.data.frame(perimeter_map)

fire_names<-as.character(perimeter_map$incidentna)
## Here we plot the shapefiles to see that they overlap
###Climate regions
plot(ExampleRaster)
plot(newfire,add=TRUE)

plot(fuel_map)
plot(newfire,add=TRUE)

Here is the processing of the Grimet data

fwi_dat_slim <- fwi_dat[,c(1,2,3,22)]
colnames(fwi_dat)

##  [1] "Year"                 "Timestep"             "EcoregionName"       
##  [4] "EcoregionIndex"       "ppt"                  "min_airtemp"         
##  [7] "max_airtemp"          "std_ppt"              "std_temp"            
## [10] "winddirection"        "windspeed"            "ndeposition"         
## [13] "co2"                  "relativehumidity"     "min_relativehumidity"
## [16] "max_relativehumidity" "specifichumidty"      "par"                 
## [19] "ozone"                "shortwave"            "temperature"         
## [22] "FWI"                  "X"

fwi_date_info <- with(fwi_dat_slim, paste(Year, Timestep))
fwi_dates <- strptime(fwi_date_info, "%Y %j") #COnverting from julian day to y-m-d
#unique(fwi_date_dat$EcoregionName)
fwi_date_dat <- cbind(as.POSIXct(fwi_dates), fwi_dat_slim[,3:4]) #attaching 
colnames(fwi_date_dat) <- c("Date", "Ecoregion", "FWI")
fwi_date_dat$Ecoregion<-gsub("eco","",fwi_date_dat$Ecoregion)
fwi_date_dat$Ecoregion<-as.numeric(as.character(fwi_date_dat$Ecoregion))

Visually filtering fires.

Looking for fires, that are not just the final burn shape.

  years <- 1992:2017
  library(raster)
  climate_day_mat <- NULL
 # One_fire <- newfire[newfire$incidentna,] #This pulls out all fire names
  
  ##Select for fires that lasted more than 1 day
  fire_names_manydays <- unique(fire_names[duplicated(fire_names)])
  fire_names_manydays <- fire_names_manydays[!is.na(fire_names_manydays)] ##Remo
  for (i in 1:length(fire_names_manydays)){
    
    #fire_select <- perimeter_map[perimeter_map$fire_name == fire_names_manydays[i],] #selecting an individual fire_ doesn't work 7/7 #
    #print(fire_names_manydays[i])
    fire_select <- subset(perimeter_map, perimeter_map$incidentna ==fire_names_manydays[i])#selecting an individual fire
    fire_days <- as.character(sort(unique(fire_select$perimeterd)))
  
   #print(max)
   #print(max[4]) 
   #if(is.na(fire_days)){fire_days<-as.character(sort(unique(fire_select$date_)))}
    #sorting fire days into chronological order
   # print(fire_days)
    ##Add a short escape in case fire days <= 1
    if(length(fire_days) < 2) next
    #par(mfrow=c(3,3))
    for(j in 1:(length(fire_days))){
    max<-raster::bbox( subset(fire_select, fire_select$perimeterd == max(fire_days)))
    fire_day_select <- subset(fire_select, fire_select$perimeterd == fire_days[j])
    if(j==1){
       plot(fire_day_select,axes=T,main=paste0(fire_day_select$incidentna,"   Day number ",j),xlim=c(max[1],max[3]),ylim=c(max[2],max[4]),col=adjustcolor('red',alpha.f = .01))
    }else(
    plot(fire_day_select,axes=T,add=T,main=paste0(fire_day_select$incidentna,"   Day number ",j),xlim=c(max[1],max[3]),ylim=c(max[2],max[4]),col=adjustcolor('red',alpha.f = .01)))
    #selecting the first fire perim from that date, in case there are multiples
    }
  }

We filtered the fires down to the final fires, based on having more than one day, and visual progress that appears to be more than mearly the same shape drawn at differnt times

SoundInc<-c("BOTELER","Halls Top","SILVER MINE",'Chestnut Knob',"Rough Ridge","DICK'S CREEK",'CLIFFSIDE','Rattlesnake Hollow','MAPLE SPRINGS','Rock Mountain',"OLD ROUGHY",
"KNOB","Quarry Creek","DOBSON 3","CLEAR CREEK","TELLICO","TIMBER RIDGE","East Miller Cove","FEREBEE MEMORIAL","HORTON","CAMP BRANCH","DOBSON KNOB","Wolf Creek","HAPPY VALLEY RIDGE","WHITE CREEK",
"Pinnacle Mountain","PARTY ROCK")

Firesinclude<-perimeter_map$incidentna[perimeter_map$incidentna %in% SoundInc]

climate_day_mat <- NULL

wsv_dat_df<-as.data.frame(wsv_dat)
##Select for fires that lasted more than 1 day
fire_names_manydays <-fire_names[!duplicated(fire_names)]
fire_names_manydays <- fire_names_manydays[!is.na(fire_names_manydays)] ##Removing potential NAs
fire_names_manydays<- fire_names_manydays[fire_names_manydays %in% SoundInc]

### For each fire for each day. 
for (i in 1:length(fire_names_manydays)){
 # print(fire_names_manydays[i])

  fire_select <- subset(perimeter_map, perimeter_map$incidentna ==fire_names_manydays[i])#selecting an individual fire
  fire_days <- as.character(sort(unique(fire_select$TrueDate)))
  ### Check that there are two days worth of data to look at. 
  if(length(fire_days) < 2) next
  storedcells<-NULL
  dftest<-as.data.frame(fire_select)
   ### This subloop looks through each set of polygons (t, t+1) and finds the last day in which expansion occurs
  ### this will indicate when to tell th algorithim to consider the remaining cells failed spread. 
  
  for(j in 1:(length(fire_days)-1)){
    ## Day t shape
    fire_day_select <- subset(fire_select, fire_select$TrueDate == fire_days[j])# selecting the first day of the fire
    fire_day_select<-fire_day_select[order(fire_day_select$acres),]
    fire_day_select <- fire_day_select[1,] #selecting the first fire perim from that date, in case there are multiples
    ## Day t+1 shape
    fire_day_select_1 <- subset(fire_select, fire_select$TrueDate== fire_days[j+1])# getting day 2 fire perimeter
    fire_day_select_1<-fire_day_select_1[order(fire_day_select_1$acres),]
    fire_day_select_1 <- fire_day_select_1[1,] #selecting the first fire perim from that date, in case there are multiples
    area_expansion <-(area(fire_day_select_1)/4046.86)-(area(fire_day_select)/4046.86)
    if(area_expansion < 30.00){
      next()
      }
    lastday<-j
 
  }
 
  for(j in 1:(length(fire_days)-1)){
 
    Failedvalues<-NA
    Successcells<-NA
    EndFailedvalues<-NA
    #print(paste0("day  ",j))
    ## day t shape
    fire_day_select <- subset(fire_select, fire_select$TrueDate == fire_days[j])# selecting the first day of the fire
    fire_day_select<-fire_day_select[order(fire_day_select$acres),]
    fire_day_select <- fire_day_select[1,] #selecting the first fire perim from that date, in case there are multiples
    ## day t+1 shape
    fire_day_select_1 <- subset(fire_select, fire_select$TrueDate== fire_days[j+1])#getting day 2 fire perimeter
    fire_day_select_1<-fire_day_select_1[order(fire_day_select_1$acres),]
    fire_day_select_1 <- fire_day_select_1[1,] #selecting the first fire perim from that date, in case there are multiples
    ## calculate area of expansion, Area m2 to acres.
    area_expansion <-(area(fire_day_select_1)/4046.86)-(area(fire_day_select)/4046.86)
    ### In this case the spread would be less than one 250m cell in our simulation.
    if(area_expansion < 30.00){
      #print("next")
      next()
      }
    
    ### Here we associate the wind and fire weather index, and reclassify the 
    ### raster to associate with spread/not spread cells
    wind_today<-wsv_dat_df %>%
      subset(date==gsub("/","-",paste0(fire_days[j],'T00:00:00Z')))
    WindMap<-climate_stack$X11_Ecoregions
    WindTrans=as.matrix(data.frame(was=c(1,2,3,4,5,6,8,9,10,11),
               is=as.numeric(t(wind_today[2:11]))))
    Wind_Reclas<-reclassify(WindMap,WindTrans)
    
    fire_today<-fwi_date_dat%>%
      subset(Date ==gsub("/","-",paste0(fire_days[j],'T00:00:00Z')))%>%
      arrange(Ecoregion)
    FireMap<-climate_stack$X11_Ecoregions
    FireTrans=as.matrix(data.frame(was=c(1,2,3,4,5,6,8,9,10,11),
               is=fire_today$FWI))
    FireWeather<-reclassify(FireMap,FireTrans)  
    End_stack<-stack(climate_stack,FireWeather)
    End_stack<-stack(End_stack,Wind_Reclas)
    
    #Creating vector of fire dates. It doesnt work when I bind them below for some reason 
    date_char <- as.POSIXct(as.character(fire_day_select$TrueDate)) 
    
    ##Extracting climate and fuels variables from raster stack.
    #This will compile a data frame of every cell that fire perimeters touch and whether 
    day1 <- extract(End_stack, fire_day_select, cellnumbers = TRUE, df=TRUE, method="simple")
    
    day2 <- extract(End_stack, fire_day_select_1, cellnumbers = TRUE, df=TRUE, method="simple")
    
    # First find day not in stored, then find adajecnts, remove adjecnets in stored. 
    
    day1cells<-day1$cell[!day1$cell %in% storedcells]
    storedcells<-c(storedcells,day1$cell)
    if( length(day1cells) >0){ 
    adjacentcells1<-adjacent(End_stack,day1cells,directions=4, pairs=TRUE,id=TRUE)
    adjacentcells1<-unique(adjacentcells1[,3])
    adjacentcells1<-adjacentcells1[!adjacentcells1 %in% storedcells]
    }else(adjacentcells1=NA)
    newspread<-day2$cell[!day2$cell %in% storedcells]
    ### Calculate failed spread, areas where spred could occur, but are have not previously been tried. 
    Failedvalues<-data.frame()
    if( length(day1cells) >0){ 
    FailedCells=adjacentcells1[!adjacentcells1 %in% newspread]
    Failedvalues<-as.data.frame(extract(End_stack,adjacentcells1))
    }
    if (nrow(Failedvalues)>0){ Failedvalues$spread_success<-0}
    Successcells<-as.data.frame(extract(End_stack,newspread))
    if (nrow(Successcells)==0){next}
    ### Successful spread gets one 
    Successcells$spread_success<-1
    ### Here if it is the last day of spread, calculate one more days possible spread,
    ### all of which is counted as failed. 
    EndFailedvalues<-data.frame()
    if(j==(lastday)){
    ### Creating the same wind and fire weather index as above. 
    wind_today<-wsv_dat_df %>%
      subset(date==gsub("/","-",paste0(fire_days[j],'T00:00:00Z')))
    WindMap<-climate_stack$X11_Ecoregions
    WindTrans=as.matrix(data.frame(was=c(1,2,3,4,5,6,8,9,10,11),
               is=as.numeric(t(wind_today[2:11]))))
    Wind_Reclas<-reclassify(WindMap,WindTrans)
    
    fire_today<-fwi_date_dat%>%
      subset(Date ==gsub("/","-",paste0(fire_days[j],'T00:00:00Z')))%>%
      arrange(Ecoregion)
    FireMap<-climate_stack$X11_Ecoregions
    FireTrans=as.matrix(data.frame(was=c(1,2,3,4,5,6,8,9,10,11),
               is=fire_today$FWI))
    FireWeather<-reclassify(FireMap,FireTrans)  
    End_stack<-stack(climate_stack,FireWeather)
    End_stack<-stack(End_stack,Wind_Reclas)
    EndCells<-adjacent(End_stack,day2$cell,directions=4, pairs=TRUE,id=TRUE)
    EndCells<-unique(EndCells[,3])
    EndCells<-EndCells[!EndCells %in% storedcells]
    EndFailedvalues<-as.data.frame(extract(End_stack,EndCells))
    EndFailedvalues$spread_success<-0
    }
    ### Adding cells to one dataframe.
    if(nrow(Failedvalues)>0 & nrow(Successcells)>0){
    Dfout<-rbind(Failedvalues,Successcells)}else{Dfout<-Successcells}
    if(nrow(EndFailedvalues)>0){Dfout<-rbind(Dfout,EndFailedvalues) }
   # print(table(Dfout$spread_success))
    Dfout$Area_expansion<-area_expansion
    climate_day_df <- cbind(fire_names_manydays[j],
                                 as.character(date_char),Dfout) #Putting everything into a big dataframe with all climate variables
    climate_day_mat <- rbind.data.frame(climate_day_mat, climate_day_df) #binding individual days to everyother day
  }
  
}
## Save for future analysis
write.csv(climate_day_mat,"Spread_3_3.csv")

climate_day_mat<-read.csv("C:/Users/zacha/Desktop/SCRPPLE_10_26/Spread_3_3.csv")
climate_day_total<-climate_day_mat[-1]
table(climate_day_total$spread_success)

## 
##    0    1 
## 6846 7368

#climate_day_total<-climate_day_mat[climate_day_mat$Area_expansion>,]


climate_day_total<-climate_day_total[climate_day_total$X11_Ecoregions.1!=0,]
#climate_day_total<-climate_day_total[climate_day_total$area_expansion>=0,]


climate_day_complete <- climate_day_total[complete.cases(climate_day_total[2:10]),]
table(climate_day_complete$spread_success)

## 
##    0    1 
## 6846 7368

#climate_day_complete<-climate_day_complete[,-1]

##Attaching a unique ID to each row in case we need it later
climate_day_complete <- cbind(1:nrow(climate_day_complete), climate_day_complete)

fix<-function(x){
return(as.numeric(as.character(x)))  
}

##Renaming columns
#climate_day_complete<-climate_day_complete[,-1]
colnames(climate_day_complete) <- c("ID","FireName",  "date", "wind_region", "fuel_number", "uphill_azi", "slope",
                               "FWI","WSPD" , "spread","expansion")
climate_day_complete$expansion<-climate_day_complete$expansion*0.404686 ## To hectares
climate_day_complete$fuel_number<-climate_day_complete$fuel_number/3000
climate_day_complete$fuel_number[climate_day_complete$fuel_number >1.0]<-1.0



U_b <- 5 #mystery value of 'combustion bouyancy'. This changes based on fire severity. Check out GItHub and how its coded in SCRAPPLE itself

# relative.wd <-  wind.direction - upslope.azi
relative_wd <- as.numeric(climate_day_complete$WSPD) - as.numeric(climate_day_complete$uphill_azi)
climate_day_complete$effective_wsv <- U_b * ((fix(climate_day_complete$WSPD)/U_b) ^ 2 + 2*(fix(climate_day_complete$WSPD)/U_b) *  
                sin(fix(climate_day_complete$slope)) * cos(relative_wd) + (sin(fix(climate_day_complete$slope))^2)^0.5)


hexbinplot(fix(climate_day_complete$spread)~fix(climate_day_complete$FWI),xlab="FWI",ylab="Spread",xbins=50,aspect=1,type="r")

hexbinplot(fix(climate_day_complete$spread)~fix(climate_day_complete$effective_wsv),xlab="Wind Speed",ylab="Spread",xbins=50,aspect=1,type="r")

hexbinplot(fix(climate_day_complete$spread)~fix(climate_day_complete$fuel_number),xlab="Fuel number",ylab="Spread",xbins=50,aspect=1,type="r")

hexbinplot(fix(climate_day_complete$expansion)~fix(climate_day_complete$FWI),xlab="FWI",ylab="Spread",xbins=50,aspect=1,type="r")

hexbinplot(fix(climate_day_complete$expansion)~fix(climate_day_complete$effective_wsv),xlab="Wind Speed",ylab="Spread",xbins=50,aspect=1,type="r")

spread_vars_short<-climate_day_complete
table(spread_vars_short$spread)

## 
##    0    1 
## 6846 7368

Full_logit <- glm(spread ~fix(FWI), 
                  data = spread_vars_short, family = "binomial")
summary(Full_logit)

## 
## Call:
## glm(formula = spread ~ fix(FWI), family = "binomial", data = spread_vars_short)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.6851  -1.1823   0.9176   1.1538   1.4485  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.892329   0.071383  -12.50   <2e-16 ***
## fix(FWI)     0.039216   0.002825   13.88   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 19686  on 14213  degrees of freedom
## Residual deviance: 19486  on 14212  degrees of freedom
## AIC: 19490
## 
## Number of Fisher Scoring iterations: 4

AIC(Full_logit)

## [1] 19490.06

FWI<-seq(0,max(spread_vars_short$FWI),.5)
xB<-exp((-0.892329)+ 0.039216*(FWI))
binomial2<-xB/(1+xB)

xB<-exp(( -0.892329  )+ 0.039216*(FWI))
binomial2min<-xB/(1+xB)

xB<-exp(( -0.892329  )+  0.039216*(FWI))
binomial2mean<-xB/(1+xB)

plot(FWI,binomial2,xlab="Fire Weather Index",ylab="Spread Probability",ylim=c(0,1),col="red")
points(FWI,binomial2mean)
points(FWI,binomial2min,col="blue")

spread_vars_short<-climate_day_complete
#table(spread_vars_short$spread)
Full_logit <- glm(spread ~fix(FWI)+fix(fuel_number), 
                  data = spread_vars_short, family = "binomial")
summary(Full_logit)

## 
## Call:
## glm(formula = spread ~ fix(FWI) + fix(fuel_number), family = "binomial", 
##     data = spread_vars_short)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.7373  -1.1810   0.9002   1.1367   1.5321  
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      -1.259443   0.087121 -14.456  < 2e-16 ***
## fix(FWI)          0.036959   0.002849  12.974  < 2e-16 ***
## fix(fuel_number)  0.694115   0.092975   7.466 8.29e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 19686  on 14213  degrees of freedom
## Residual deviance: 19430  on 14211  degrees of freedom
## AIC: 19436
## 
## Number of Fisher Scoring iterations: 4

#hist(spread_vars_short$finefuels)
AIC(Full_logit)

## [1] 19436.07

FWI<-seq(0,max(spread_vars_short$FWI),.5)
xB<-exp(( -1.259443 )+ 0.036959*(FWI)+0.694115*max(spread_vars_short$fuel_number))
binomial2<-xB/(1+xB)

xB<-exp(( -1.259443 )+ 0.036959*(FWI)+0.694115*min(spread_vars_short$fuel_number))
binomial2min<-xB/(1+xB)

xB<-exp(( -1.259443 )+ 0.036959*(FWI)+0.694115*mean(spread_vars_short$fuel_number))
binomial2mean<-xB/(1+xB)

plot(FWI,binomial2,xlab="Fire Weather Index",ylab="Spread Probability",ylim=c(0,1),col="red",type="l",lwd=3.0,cex=1.2,cex.axis=1.2,cex.lab=1.2,
     main="Probablity of Spread")
lines(FWI,binomial2mean,lwd=3.0)
lines(FWI,binomial2min,col="blue",lwd=3.0)
legend(0,1.0,legend=c("High Fine Fuels","Median Fine Fuels"," Low Fine Fuels"),
       lty=c(1,1,1),lwd=c(3,3,3),col=c("red","black","blue"))

Max spread.

Take the max spread per FWI and WS fit to that.

check<-aggregate(x=list(Ws=climate_day_complete$effective_wsv,FWI=as.numeric(climate_day_complete$FWI)),
          by=list(date=climate_day_complete$dat,expansion=as.numeric(climate_day_complete$expansion)),FUN=mean)

check$RWS<-round(check$Ws)
check$RFWI<-round(check$FWI)



mean(check$expansion)

## [1] 409.1673

plot(fix(check$expansion)~fix(check$FWI),xlab="FWI",pch=19,col=adjustcolor("red",alpha.f = .3),ylab="Spread",cex.lab=1.2,cex.axis=1.2)

plot(fix(check$expansion)~fix(check$Ws),xlab="Wind Speed",ylab="Spread",pch=19,col=adjustcolor("red",alpha.f = .3),cex.lab=1.2,cex.axis=1.2)

Expansionfull<-lm(fix(expansion) ~FWI , data =check)
#?glm()
summary(Expansionfull)

## 
## Call:
## lm(formula = fix(expansion) ~ FWI, data = check)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -480.4 -301.7 -146.2  146.6 3248.5 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  282.877    206.852   1.368    0.174
## FWI            5.059      8.086   0.626    0.533
## 
## Residual standard error: 496.6 on 119 degrees of freedom
## Multiple R-squared:  0.003278,   Adjusted R-squared:  -0.005098 
## F-statistic: 0.3914 on 1 and 119 DF,  p-value: 0.5328

FwI<-seq(0,max(spread_vars_short$FWI),.5)
FWI_Spread<- 282.877+(5.059*FwI)
plot(FwI,FWI_Spread,xlab="FWI",ylab="Spread per day")

None of the variables were predcitive to maximum fire size: We ended using the mean of the fire spread

mean(check$expansion)

## [1] 409.1673